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PREFACE 

A  basic  problem  in  radiative  transfer  is  the  estimation  of  the 
physical  parameters  of  a  scattering  and  absorbing  medium  based  on 
measurements  of  the  diffusely  reflected  light.  In  this  series  of 
papers  we  shall  show  how  such  tasks  may  be  viewed  .as  nonlinear 
boundary-value  problems  which  may  be  solved  computationally  using 
the  technique  of  quasilinearization. 

In  this  Memorandum  we  consider  a  stratified  medium  consisting 
of  two  layers.  Our  aim  is  to  determine  the  optical  depth  and  the 
albedo  for  single  scattering  of  each  layer  based  on  measurements  of 
the  angular  dependence  of  the  light  diffusely  reflected  from  the  slab. 
The  results  of  some  numerical  experiments  are  presented,  and  FORTRAN 
IV  listings  are  provided. 

The  results  will  be  of  interest  to  meteorologists,  experimental 
physicists,  and  control  engineers. 
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SUMMARY 

Confider  a  slab  which  consists  of  two  layers  which  absorb  and 
scatter  light.  First  the  basic  differential-integral  equation  for 
the  intensity  of  the  diffusely  reflected  light  is  given,  the  source 
being  uniform  parallel  rays.  Then  we  show  how  to  determine  the 
nature  of  each  layer  using  measurements  of  the  light  diffusely  re¬ 
flected  from  the  slab.  We  view  this  as  a  nonlinear  boundary- value 
problem  and  show  that  it  may  be  resolved  computationally  using  the 
technique  of  quasilinearization.  Numerical  examples  and  FORTRAN 
programs  are  provided. 
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I.  INTRODUCTION 

Problems  of  radiative  transfer  in  planetary  and  stellar  atmos- 

(12  3) 

pheres  have  been  extensively  treated  by  many  investigators. '  »  * 

Much  of  this  work  has  dealt  wi  h  what  we  may  call  direct  problems, 
i.e.,  the  determination  of  the  intensities  of  radiation  produced  by 
certain  light  incident  on  a  medium  with  known  scattering  and  absorbing 
properties.  Recently,  though,  as  a  growing  list  of  papers  attests, 
there  has  been  increasing  interest  in  inverse  problems.  In  these, 
the  aim  is  to  deduce  the  nature  of  the  medium  (and  often  the  source) 
on  the  basis  of  measurements  of  the  radiation  field  produced  by  in¬ 
cident  radiation. 

With  this  Memorandum  we  initiate  the  study  of  inverse  problens 

(7  g\ 

in  radiative  transfer  via  the  method  of  quasilinearization. ^  ^  We 

wish  to  show  that  such  problems  can  be  viewed  as  multipoint  boundary- 
value  problems  for  large  systems  of  nonlinear  ordinary  differential 
equations.  These  equations  can  then  be  resolved  numerically  using 
current  digital  computers  and  the  method  of  quasilinearization. 

We  consider  a  slab  which  consists  of  two  layers  of  isotropically 
scattering  material  of  unknown  optical  thickness  and  albedo  for  single 
scattering.  Given  the  measured  intensity  of  the  light  diffusely  re¬ 
flected  from  the  slab,  we  wish  to  determine  the  optical  thickness  of 
the  two  layers  and  their  albedos  for  single  scattering. 

Ip  subsequent  papers  we  shall  consider  more  complex  processes 
involving  ani'jotropic  scattering,  time-dependence,  and  noncoherent 


scattering. 


We  consider  an  inhomogeneous,  plane-parallel,  non-emitting  and 
isotropically  scattering  atmosphere  of  finite  optical  thickness  Tj 
whose  optical  properties  depend  only  upon  t»  the  optical  height  above 
the  bottom  (0  S  t  S  .  Let  parallel  rays  of  light  of  ne*"  flux  tt 
per  unit  area  nonnal  to  their  direction  of  propagation  be  incident 
on  the  upper  surface  in  the  direction  characterized  by  the  number 

(0  <  u  <  1),  where  u.  is  the  cosine  of  the  angle  measured  from 
o  o  o 

the  normal  to  the  surface.  The  bottom  surface  is  a  completely  ab¬ 
sorbing  barrier,  so  that  no  light  is  reflected  from  it.  This  is  not 
an  essential  requirement.  See  Fig.  1. 


Fig.  l-~ Incident  and  reflected  rays  for  a  slab  of  thickness  t, 
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The  intensity  of  the  diffusely  reflected  light  in  the  direction 
cos  ^  p.  is  r(p.,  Tj^).  We  define  a  related  function  R(ii,  '’’P  * 
syimetric  in  and  by  writing 


M-qJ  t^) 


R(^if  T^) 

r; 


(1) 


The  function  R  satisfies 


the  integro-dif ferential  equation 


(1,9) 


SR(^i,^ip, 


.  ri 


^  )  Mm., 


(2) 


where  X(Tj^)  is  the  albedo  for  single  scattering,  and  R  is  subject  to 
the  condition 

-  0.  (3) 


This  equation  is  obtained  by  means  of  the  theory  of  invariant  imbed** 
ding.  A  discrete  version  of  Eqs.  (2)  and  (3)  is  obtained  by  replacing 
the  integrals  by  finite  sums  using  a  Gaussian  quadrature  formula  of 
order  N.  Let  us  introduce  new  functions 


R^jCTi)  “  i,j  »  1,2,. ..,N, 


(4) 


where  1®  roots  of  the  shifted  Legendre  polynomial 

* 

Pjj(p)  =  Pjj(l-2^i),  The  values  of  and  the  corresponding  Christoffel 
weights  are  tabulated  in  Ref.  10  for  N  -  3, 4,..., 15.  The  ftmctions 
Rij(Ti)  satisfy  the  system  of  ordinary  nonlinear  differential  equation 
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X(Tj)  [l-x  I  I  J  ]  [l+  -1  I  ^  ]  . 


(5) 


with  the  initial  condltiona 

Rij(O)  -  0  ,  (6) 

for  i  -  1,2,. ..,N;  J  -  1,2,. ..,N. 


5 


III.  AN  INVERSE  PROBLEM 


Let  us  now  consider  a  medium  composed  n£  two  layers,  with  albedo 
for  single  scattering  in  the  lower  layer  and  albedo  X2  upper 

layer  (Xj^  ^  X2^  •  total  optical  thickness  of  this  slab  is  T  and 

the  thickness  of  the  lower  layer  is  c.  See  Fig.  2.  We  wish  to  de¬ 
termine  c,  T,  and  the  parameters  Xj^  and  X2- 


Fig.  2— A  stratified  medium 


Let  us  assume  that  we  have  obtained  N^  (noisy)  measurements  of 
the  diffusely  reflected  light,  whe*:e  b^^  is  the  intensity  of  the 

light  diffusely  reflected  in  the  direction  caused  by  Incident 
parallel  rays  of  net  flux  rr  in  the  direction  The  constants  which 

characterize  this  medium,  Xj^»X2»®>  and  T  are  to  be  determined  so  that 


the  theoretical  diffuse  reflection  pattern  produced  by  using  the 
estimated  values  in  the  differential  Eqs.  (5)  will  agree  as  closely 
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as  possible^  In  the  least  squares  sense,  with  the  observed  field 


(bijJ,  i.e. ,  we  wish  to  minimize  the  expression 

f  i  i- 

M  j»l 


(7) 


Let  us  assume  that  the  albedo  function  has  the  form 


X(t)  ■  a  +  b  tanh  10  (t  -  c) 


(8) 


so  that 


^  a  -  b  in  Layer  1, 

X2  —  A  +  b  in  Layer  2, 


(9) 


and  c  is  the  position  of  the  boundary  between  Layer  1  and  Layer  2. 

See  Fig.  3. 

The  "observations"  {8^^)  are  produced  computationally  with  the 
use  of  Eq.  (8)  for  the  albedo  function,  where  we  set 

a  «  0.5,  b  «  0.1,  c  -  0.5,  (10) 

and  with  the  use  c£  the  differential  Eqs.  (5),  integrating  out  to  a 
thickness 

T  =  1.0  . 


The  observations  which  this  produces  are  given  in  Table  1. 


T  =  1  .0 


Fig.  3— The  albedo  function 


Table  1 

THE  MEASUREMIOTS  fb.  J 


■n 

■Hi 

2 

3 

4 

5 

6 

7 

i  =  1 

0.079914 

0.028164 

0.014304 

0.009104 

0.006707 

0.005515 

0.004970 

CM 

0.143038 

0.091522 

0.058437 

0.040826 

0.031405 

O.C26378 

0.023989 

3 

0.167000 

0.134331 

0.099653 

0.075106 

0.060044 

0.051445 

0.047248 

4 

0.178898 

0.157955 

0.126408 

0.099392 

0.081253 

0.070435 

0.065042 

5 

0.185284 

0.170817 

0. 142072 

0.114229 

0.094495 

0.082423 

0.076332 

6 

0.188723 

0.177733 

0.150791 

0.122665 

0.102104 

0.089349 

0.082870 

7 


0. 190354 b. 180898  0. 154995  0. 126773  0. 105829  C. 092748  0.086083 
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IV.  NllMERICAI.  SOLUTION  VIA  QUASILINEARIZATION^^*^* 

Let  us  consider  the  problem  posed  above  in  more  general  and 
flexible  terms.  An  R- dimensional  column  vector  x(t)  is  a  solution 
of  the  differential  equation 

X  -  f(x)  ,  0  S  t  ^  T  .  (11) 

Tlie  value  of  T  ."'s  assumed  known.  This  again  represents  no  essential 
restriction  as  we  shall  see.  We  wish  to  determine  the  initial  vector 

c, 

x(0)  =  c  ,  (12) 

in  such  a  manner  that  we  minimize  the  quadratic  form 

(x(T)  -  b,  x(T)  -  b)  -  Q  ,  (13) 


where  x(T)  is  the  solution  of  Eq.  (1)  for  t  ■  T.  The  first  a  com¬ 
ponents  of  c  are  required  to  be  zero  and  the  remaining  ones  are 
variables  to  be  determined. 

Our  computational  formalism  proceeds  as  follows:  First  an  initial 

approximation  to  the  desired  initial  vector  is  selected.  We  then 

proceed  inductively.  Suppose  that  we  have  obtained  a  approxima- 

k  th 

tion  to  c,  which  we  denote  by  c  ,  and  a  k  approximation  to  the 

solution  function  x(t),  which  we  denote  by  x  (t).  Note  that  the 

superscripts  refer  to  the  order  of  the  approximation  and  not  to  the 

k+1  k+1 

components  of  c  and  x.  We  obtain  the  vectors  c  and  x  (t)  in 

k+1 

this  manner.  Tlie  vector  x  (t)  is  a  solution  of  the  linear  system 
k+1 


in  X 


X  -f(x)+J(x)(x  -x), 


(14) 
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where  J(x  )  is  the  Jacobian  matrix  with  elements 


-  (Bf^/ax^)  .  (15) 

We  produce  a  particular  solution  p(t)  of  the  system  (14)  by  numerically 
integrating  the  system 

p  -  f(x*')  +  J(x*')  (p  -  x*^)  ,  (16) 

where  p(t)  Is  subject  to  the  initial  condition 


p(0)  »=  0  . 


(17) 


Then  we  produce  numerically  R-s  independent  solutions  of  the  homoge¬ 
neous  system 

h.  -  J(x'')h^  ,  (18) 

where  the  vector  h^(t)  is  subject  to  the  initial  condition 

(j*”^  component  of  h^(0))  -  6^^  ,  i  ■  s  +  1,  R  ,  (19) 


being  the  Kronecker  delta  function.  These  integrations  are  readily 
carried  out  on  the  interval  0  S  t  s  T,  since  complete  sets  of  initial 
conditions  are  given.  From  general  theory  we  know  that  the  vector 
x^^^  is  representable  in  the  form 


k+1 


(t)  =  p(t)  +  ^  nj^h^(t)  , 


(20) 


i“afl 


where  the  numbers  are  arbitrary  constants.  In  view  of  the  way  we 
have  chosen  the  initial  conditions  in  Eqs.  (17)  and  (19),  we  see  that 
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reprcKentc  the  initiel  value  of  the  component  of  the  vector 
k+1 

X  .  In  particular,  when  the  integrations  are  completed,  we  shall 

have,  in  numerical  form,  the  values  of  the  vectors  p(T),  h  (T) 

8*^  X  ^ 

We  now  wish  to  minimize  the  form 


R  R 

(p(T)  +  y  m^h^(T)  -  b,  p(T)  +  ^  m.h^(T)  -  b) 

i-^1  i«8+l 

(21) 


over  the  values  of  m^,  i  »  s+l,  8+2,  ...,  R,  where  b  is  known  experi¬ 
mentally  and  p(T),  and  h 

9+1 

putationally.  This  minimization  is  done  by  solving,  numerically, 
the  linear  alg«‘braic  system  of  equations 


(T) ,  ****  known  com- 


hf/aa^  ■  0»  f  ■  »+l»  s+2,  R  . 


(22) 


The  values  of  m 


s+1’  ®s^2’  ®R 


so  obtained  are  the  values  of  the 


k+1  k+1 

last  R-s  components  of  the  vector  c  ,  and  x  (t)  can  be  determined 
from  Eq.  (20)  using  the  values  of  *'“8+2’  stored 

values  of  p(t)  and  h^(t),  i  “  s+1,  s+2,  R.  The  entire  process 

is  then  repeated  to  obtain  the  (k+2)"*^  approximations. 

To  simplify  the  procedure  of  finding  the  m' s,  let  us  observe 
that  the  relation 

k+1 


(T) 


(23) 


can  be  written  in  the  form 
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R 

^  h.(T)in^  s-b  -  p(T)  (24) 

i=s+l 

tH 

If  we  let  H  be  the  matrix  whose  i  column  is  h^(T),  i  =  »+l ,  »f2, 

•  ••>  y  bh  column  vector  b  -  p(T)j  and  m  be  the  column  vector 
whose  i  component  is  m^,  then  the  above  relation  can  be  rewritten 

-  y  .  (25) 

where  H  is  a  matrix  with  R  rows  and  R-s  columns,  and  m  and  y  are  R-s 
dimensional  column  vectors.  According  to  the  theory  of  the  method 
of  leact  squares  the  vector  m  which  minimizes  the  form 

(Hm  -  y,  rtn  -  y)  «  G  (26) 

is 

m  -  (H^H)"^H^-  .  (27) 

where  is  the  transpose  of  H, 

A  few  general  conments  on  the  procedure  are  now  in  order.  The  se¬ 
lection  of  a  good  initial  appro':imation  is  Important,  since  in  this 
case  the  method  is  rapidly  convergent,  with  the  number  of  correct  digits 
approximately  doubling  with  each  additional  step;  if,  however,  the  ini¬ 
tial  approximation  is  too  poor,  the  method  may  be  divergent.  At  each 
step  we  must  integrate  R(R-s+l)  first  order  differential  equations  with 
given  initial  conditions  to  produce  the  R-s  homogeneous  solutions  and 
the  one  particu’ ir  solution.  In  addition  we  must  solve  a  system  cf  R-s 
linear  algebraic  equations.  This  may  be  a  sizable  computing  load,  and 
the  solving  of  the  linear  algebraic  equations  can  be  a  source  of  great 
difficulty. 
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V.  NIMERICAL  EXPERIMENTS  TO  DETE^INE  c, 
THE  THICKHESS  OF  THE  LOWER  LAYER 


In  the  first  series  of  experiments  we  wish  to  determine  only 
the  thickness  of  the  lower  layer,  assuming  that  all  of  the  other 
parameters  of  the  medium  are  known.  We  use  a  seven-point  Gaussian 
quadrature,  so  that  N  ■  7.  We  consider  c  to  be  a  function  of  for 

which  dc/dTj^  ®  0.  Following  the  method  prescribed  in  the  preceding 

2  k+1 

section,  we  obtain  IT  ■  49  linear  differential  equations  for 

k-fl 

and  one  equation  for  c  .  Thus  there  are  50  linear  differential 
equations 


IR  "  u 

1  i  J  i-1  f 


N  .  W 


^  A-1  * 


^  i-i  ^ 

X  (-  10  b  sech^  10  (t,  -  0*^))!  , 
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where 


XCc'^)  ■  •  +  b  c«nh  10(tj^  -  c^) 


Cne  may  reduce  the  number  of  differential  equations  in  the 


system  (28)  by  using  the  synnetry  property 


j  -  1,2,. ...N  . 


2  k+1  k+1 

Thus  instead  of  N  +  1  «  50  differential  equations  for  and  c  , 

we  have  (N(lffl)/2)  +  1  ■  29  differential  equations  plus  the  finite 

Eqs.  (29)  to  fill  in  the  missing  values  of  R.  However,  we  will  itill 

speak  of  the  matrix  having  7  rows  and  7  columns. 

k+1 

Now  let  the  50- dimensional  vector  x  (tj^)  represent  the  49 
k+1 

elements  R^^  (tj^)  taken  in  some  order, 


k+1.  .  „k+l.  . 

^  ■  ■'ij  '■'!> 


for  T-  1,2,. ..,49  as  i  ■  1,2,. ..,7,  and  j  «  1,2,. ..,7,  and 


k+1  .  .  k+1  .  ^ 

X50  (t^  -  c  (T^  . 


We  express  x  (t^)  as  a  sum  of  a  particular  solution  p(t)  and  a 
homogeneous  solution  h(t)  of  the  Eqs.  (28), 


X  (Tj^)  -  p(t^)  +  m  h(Tj^)  . 


We  require  that  the  multiplier  m  be  chosen  to  minimize  the  expression 
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I  [p,(X)  +  m  h^(l)  -  b J  ,  (33) 

/-I 

where  the  singly  subscripted  set  {t>^)  (I  *  1,2,..., 49)  represents 
the  doubly  subscripted  set  of  constants  "  1»2,...,7, 

j  ■  1,2,..., 7)  taken  in  the  same  order  as  the  transformation  (30). 

The  value  of  m  which  minimizes  (33)  is  that  for  which 


m 


h^(i)[b^  - 

49  2 

2  [h,(l)3 

/-I  ^ 


(34) 


as  we  see  from  a  simple  differentiation.  With  the  proper  choice  of 
initial  values  p(0)  and  h(0) ,  this  constant  m  is  the  initial  value 


m 


k+1 

c 


(0) 


and  hence  gives  directly  the  thickness  of  Layer  1. 

The  results  of  three  experiments  with  different  initial  guesses 
of  the  thickness  c  are  shown  in  Table  2.  Th  itial  approximation 
is  generated  in  each  run  by  integrating  the  nonlinear  Eqs.  (5)  with 
the  value  of  c  listed  as  approximation  zero. 


Table  2 

SUCCESSIVE  APPROXIMATIONS  OF  c 


Approximation 

Run  i 

Run  2 

Run  3 

0 

0.2 

0.8 

0.0 

1 

0.62 

0.57 

2 

0.5187 

0.5024 

No 

3 

0.500089 

0.499970 

converge  ice 

4 

0.499990 

0.499991 

True  Vnlue 

0.5 

0.5 

0.5 

15 


Even  with  60  per  cent  errors  in  the  initial  approximations,  the 
value  of  c  is  determined  to  one  part  in  fifty  thousand  or  0.002  per 
cent  in  four  iterations.  The  time  required  on  the  IBM  7044  is  1% 
minutes  per  run,  using  an  integration  step  size  of  0.01  and  the  Adams- 
Moulton  methoa.  Making  use  of  the  synmetry  of  P,  58  linear  differen¬ 
tial  equations  have  to  be  integrated  for  the  particular  and  homogeneous 
solutions  for  each  approximation.  In  Run  3  the  solution  diverges  be¬ 
cause  the  initial  guess  is  not  good  enough. 
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VI.  NUMERICAL  EXPERIMENTS  TO  DETERMINE  T. 

THE  OPTICAL  THICKNESS 

Next  we  consider  the  estimation  of  the  total  optical  thickness  T 
when  all  of  the  other  system  constants  are  known  and  49  measurements 
tb}  are  given.  The  formulation  of  the  linear  boundary- value  problem 
proceeds  as  for  the  previous  case  with  one  major  difference.  The 
terminal  boundary  for  the  interval  of  integration  is  unknown.  In 
order  that  the  interval  be  fixed,  ie  define  a  new  independent  vari¬ 
able  a, 

^  T  -  , 

with  0  s  o  S  1  as  the  range  of  integration.  The  constant  T  is  a 
solution  of  the  equation  dT/do  »  0.  Then  one  may  obtain  the  linear 
equations  which  correspond  to  (28),  with  independent  variable  a 
instead  of  Tj^. 

Three  trials  are  made  to  determine  T,  with  initial  guesses 
T°  ■  0.9,  1.5  and  0.5  respectively.  Recall  that  the  true  value  is 
T  -  1.0.  Within  four  minutes  of  computing  time,  four  iterations  are 
carried  out  per  trial,  and  the  results  obtaineo  are  correct  to  one 


part  in  100,000. 
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VII.  NUMERICAL  EXPERIMENTS  TO  DETEroilNE  and  c. 

THE  ALBEDOS  AND  OPTICAL  THICKNESS 

In  the  final  experiment,  on  the  basis  of  the  49  observations, 
we  try  to  determine  X2>  c  assuming  T  *  1.0  is  known.  Again 
we  use  N  =  7.  This  time  there  are  three  homogeneous  solutions  and 
a  particular  solution  to  compute,  each  with  N(N+-l)/2  +  3  =  31  com¬ 
ponents,  so  that  4  X  31  “  124  linear  differential  equations  are  to 
be  integrated  during  each  stage  of  the  successive  approximation 
calculations.  A  standard  (Gaussian  elimination)  matrix  inversion 
procedure  is  used  to  invert  the  3x3  matrix  of  the  linear  algebraic 
system.  The  results  in  Table  3  are  obtained  in  three  iterations 
which  consume  2  minutes  of  computing  time  on  the  IQl  7044  machine. 


Table  3 

SUCCESSIVE  APPROXIMATIONS  OF  X^,  X2»  AND  c 


Approximation 

^1 

^2 

c 

0 

0.51 

0.69 

0.4 

1 

0. 4200 

0.6052 

0.5038 

2 

0.399929 

0.599995 

0.499602 

3 

0.399938 

0.599994 

0.499878 

True  Values 

0.4 

0.6 

0.5 
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VIII.  DISCUSSION 

The  »,v.neral  conceptual  and  computational  approach  to  inverse 
problems  which  we  have  discussed  here  is  by  no  means  limited  to  the 
simple  physical  model  considered.  At  the  expense  of  additional  com¬ 
puting  time  we  may  treat  atmospheres  having  anisotropic  scattering, 
time- dependent  sources,  and  so  on,  v»e  may  also  wish  to  consider 
wave  rather  than  particle  processes.  Applications  to  orbit  deter- 
inin«ition, system  identification, cardiology and  other 
areas^^^^  have  been  made.  Of  particular  importance  is  the  question 
of  the  effect  of  errors  in  the  observations  on  the  accuracy  of  the 
estimates  of  the  parameters.  In  this  connection  it  appears  that 
the  use  of  the  min-max  criterion  rather  than  that  of  least  squares 
may  be  efficacious.  These  and  other  matter'  will  be  discussed  in 
later  papers. 
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APPENDIX 


This  Appendix  lists  the  FORTRAN  IV  programs  which  were  written 
for  (I)  The  Estimation  of  c,  (II)  The  Estimation  of  T,  and  (III)  The 
Estimation  of  and  c.  Tne  major  difference  in  the  three 

programs  lies  in  the  subroutine  DAUX,  in  which  the  right-hand  sides 
of  the  differential  equations  being  integrated  must  be  evaluated  in 
the  same  order  as  the  dependent  variables  ther selves  and  both  sets 
stored  in  the  array  called  "T".  T(2)  contains  the  current  value  of 

the  independent  variable,  T(3)  contains  the  integration  step  length, 
the  next  M  locations  from  T(4)  through  T(M+3)  contain  the  M  variables 
which  are  being  integrated,  followed  by  the  M  derivatives  in  locations 
T(M+A)  through  T(2M+3).  The  integration  routine  used  (RAND  7044 
Library  Routine  Ntnnbei  W031)  is  called  upon  with  the  statements 
CALL  INTS  (T,M, 2, 0,0, 0,0, 0,0)  and  CALL  INIM,  and  these  in  turn  call 
upon  subroutine  DAUX. 

In  program  (III),  the  matrix  inversion  and  linear  algebraic 
equation  solving  subroutine  MATINV  (RAND  Library  Routine  Number 
W019)  is  used.  The  matrix  of  coefficients  in  the  program  is  EliAT, 
the  number  of  unknowns  is  three,  and  the  right-hand  vector  is  FVEC. 


< 

< 

< 


( 


r»r»  rt  rt  r»r»r*r» 


rinr»r»  ^  orno 


PROGRAM  FOR  TIE  ESTIMATION  OF  c,  THE  THICKNESS  OF  THE  LOtfER  LAYER 

20 

tIBFTC  RTINV 

COMMON  N.RT(7!»WT(7!»WR(7!»ARI7»7).NPPN7,MlMAX»KMAX.0ELTAtXTAU 

1  2ERLAM.XLAM(2) ,02(7.7)»R2(7»7J. I F LAG . R ( 2 8 . 1 C 1 5 . T ( 149 1 ) » S I G . 

2  P(28»101)»H(28.3»101! tPIAUtPLAM ( 2  J  *HTAUJ  3  >  tHLAM( 2*3 )  .P2 ( 7.7 ) 

3  H2(7.7.3)  .CONSTO)  .NFQ 

PHASF  I 

1  RFAOIOOO.N 
PRiNT899 
PRINT900.N 

RtADlOOl , ( RT ( I ) . 1=1 ,N) 

PRINT901  .  rRT(  I  )■»  1  =  1  .N) 

READlOOl . ( WT(  1  )  .  1  =  1  .N) 
pr:nt90i. ( wt ( 1 ) . 1=1 .n) 

DO  2  1=1 .N 

WR(  I  ) =WT( I ) /RT( I ) 

DO  2  J=1.N 

2  AR(1.J)=  l.O/RKI)  +  1*0/RT(J) 

899  FORMAT ( 1H146X36HRA01AT IVE  TRANSFER  -  INVERSE  PROBLEM  /  ) 

1000  F0RMAT(6I12) 

900  FORMAT(6I20) 

1001  F0RMAT(6E12.8) 

901  FORMAT{6E20.8) 

REAOICOO.NPRNT.MIMAX.KMAX 
PRINT900.NPk’NT,M1MAX,KMAX 
READlOOl. delta 
PRINT901, DELTA 

READlOOl ♦ XT AU.2ERL AM. X LAM ( 1 ) .XLAM(2) 

PRINT902 

PR1NT903,XTAU.2ERLAM.XLAM( 1 ) .XLAM(2 ) 

902  FORMAT! lH123HPHASe  I  -  TRUE  SOLUTION  /) 

903  format (IHO/ 

1  IXIIHTHICKNESS  =.  F10.4  / 

2  IXllHALBEDOl X)  «»  20HA  +  B*T ANH ( 1 0* ( X-C ) )  // 

3  1X3HA  *,  E16.8.  10X3HB  «.  E16.8,  10X3HC  ».  E16.0  //) 

call  nonlin 

00  3  1=1,N 
DO  3  J*1.N 

3  02(1 »J)«R2( I tJ) 


\ 

PHASE  I  I 

4  READlOOl .XTAU.ZERLAM.XLAM I  1 ) »XLAM( 2 ) 
K  =  C 

PRINT904.< 

PR  I NT903, XTAU.ZERLAM.XLAM! 1 ) .XLAM12  ) 
C 

CALL  NONLIN 
C 

904  F0RMAT!1H1  1 3HAPPR0X I MA T I  ON .  13/  ) 


C 

C 


OUASILINEARIZATION  ITERATIONS 


n  o 


C 


21 


DO  5  Kl=l iKMAX 
PRINT904iKl 
CALL  PANDH 

call  linear 

5  CONTINUE 


READlOOO, IGO 
GO  TO  I  1 »4)  1  IGO 
END 

5IBFTC  DAUX 

SUBROUTINE  DAUX 

DIMENSION  V2(7»7)  tXO)  »F{7)  *0(7) 

COMMON  NfRT<7) ♦WT<7) »WR<7) .AR<7.7) . NPRNT .MIMAX . <MAX » DEL TA f X TAU 

1  ZERLAMf XL.AM( 2 )»B2(7,7).R2<7.7)» I  FLAG »R I 28»101)»T(1491)iSIGf 

2  P(28fl0l)fH<28.3»101) . P T AU . PLAM ( 2 ) .HTAU<  3 ) ♦HLAM(2  »3J  tP2  <  7»7) 

3  H2(7,7,3) »C0NST<3) tNEQ 
GO  TO  ( 1»2 )  ♦  IFLAG 

C 

CNONLINEAR 

C 

1  L  =  3 

DO  4  I=1,N 
DO  4  J=1,I 
L  =  L  +  1 

4  V2(I»J)=T(L1 
DO  5  1=1. N 
DO  5  J=I.N 

5  V2 ( I ♦J|=V2< J»  I  I 
L  =  L  +  1 

VLAM2=T(L) 

SIG=T<2) 

Y=XTAU*SIG 
X ( 1 ) sZERLAM 
X<2)=XLAM( 1) 

X(3)=VLAM2 

CALL  AL0EDO( Y.X.Z) 

ZLAMDA=2 

C 

DO  6  I  =1 »N 
F  <1 ) =0.0 
DO  7  K=1,N 

7  FU  )=Fn  )  +  WR<K)*‘v2(  I  »K  ) 

6  F< 1 )=0.5*F(  I  )  +  1.0 
C 

DO  8  1  =  1  .N 
DO  8  J=1.I 
L  =  L  +  1 

0R=-ARI I ♦ J)*V2( I .J)  +  ZLAMOA*F( I )*F<J) 

8  T(L)=DR 

DO  9  1=1.1 
L  =  L+1 

9  T<L)*0.0 


RETURN 
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C 

C 

CLINEAR 

C 

2  SIG*T{2) 

Y  =  XTAU*.SIG 
X ( 1  )=ZERLAM 
X{2)=XLAM{  1  ) 

X(3)*XLAM{2) 

CALL  AL8EDC( Y.X»Z) 

ZLAMDA*Z 

/* 

V. 

DO  16  1=1 tN 

r{  I  )=o.o 

DO  17  K=1»N 

17  F{  I  )*F{ I )  +  WR{<)*R2(  I »K  ) 

16  F{  I  )*0.5»F{ I )  +  l.O 
C 

CP»S 

C 

L  =  3 

DO  14  1 « 1 ♦ N 
DO  14  Jaltl 
L  =  L  +  1 

14  V2{I*J)»T(L) 

DO  15  1*1 »N 
DO  15  J=I ,N 

15  V2( I  .J)=V2IJ* I ) 

L  =  L+1 

VLAM2=T{L) 

C 

DO  10  1=1. N 
Gl  I  )=0.C 
DO  10  <»1.N 

10  G(I)=G(I)  +  {  V2<  I  »K)-R2{  I  )  )*WR()C) 

ARG=10.C*{Y-XlAM(2)  ) 

XTANX  =  -10.0*XLAM{ 1 )*( 1 , 0- <  T ANH { ARG )  )**2) 
M=3+NEQ 
DO  1 2  I  =  1  .  N 

no  1?  j=i.i 

FIJ»F{ I )*F( J) 

CAPF  =  -AR(  I  .  J)  *R2  {  I  »  J  )  +  ZLAN‘OA*FI^ 

Tl'x  -AR(  I  .J)*{V2(  I  .J)-R2n  .J)  ) 

T2«  0. 5*ZLAMDA*{F{ I ) *G ( J ) +F ( J ) *G ( I ) ) 

T3=  CAPF 

T4= ( V^AM2-XlAM{ 2 ) )»XTANX*FIJ 
M  =  M+1 

12  T  (M)  =Tl>t>T2  +  T3+T4 
DO  19  1=1.1 
M  =  M+1 

19  T(.y)=0.0 

c 

CH'S 


c 


DO  100  <=1,1 
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00  24  1*1, N 
00  24  J=1,I 
L  =  L  +  1 

V2 ( I , J)=7(l) 

DO  25  1=1. N 
DO  25  J=i.^^ 

25  V2(I,j).v2(j,,, 

L  =  Ltl 

VLAM2  =  T(L  > 


c 
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DO  20  I  =  i,»^ 

G( I ) =0.0 
DO  20  J.i  ,N 

G'I)=G(n  +  V2(I.J)*WR(J, 


DO  22  1  =  1, M 
DO  2?  J=i.r 
f'ijsF  ( n#F  ( j) 

I»J)*V2(I,J) 

Tj'o.c 

T4k VLAM2»XTANX#F IJ 

M*M+1 

^  22  T(M).ti+T2+T3+T4 


29 

100 


DO  29  r-1,1 
M  =  M+1 

T ( M ) =0 .0 
CONTINUE 
return 

END 

SIBFTC  NONLIN 

subroutine  nonlin 

COMMON  N»RT(7)  WT/- 

1  ^^•^>-AM»XLAM(2  )  ,  Bp  !  ^  ^  ’  *^RRNT  .MIMAX  ,  <WA  y  riFi  TA 

'"”IAL  APPROX. 

T(2)=Q.0 

T(3)=DELTA 

M=1 

L1=0 


L'>  =  3 

DO  1  I=i,N 
DO  1  J=i,i 
Ll=Li+l 
L3*L3+1 

R2 ( I » J )»o.O 

R(L1»M)=R2( r .J, 
1  T{L3) =R2( I  ,J) 

L3=L3+1 


U  KJ  KJ  KJ 
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2  T{L3) =XLAM{2) 

C 

NPO» { N*<  N+1 )  ) /2  ♦  1 

CALL  INT5{ T»NEQ»2»0»0*0»0*0»0) 

C 

SIG*T(2) 
call  OUTPU’i 
C 

DO  5  M1.1,M1MAX 
DO  4  M2=1»NPRNT 
call  INTM 

M  =  M+1 
LI  =0 
L3»3 

DO  3  I=1»N 
DO  3  J=1»I 
I  1 =L1  +  1 
L3=L3+1 
R2{ I ♦J)=T{L3) 
b  R(L1 tMjsRZJ I » J) 

4  SIG=T{2) 

5  CALL  OUTPUT 
C 

RETURN 

END 

SIBFTC  PANOH 

SUBROUTINE  PANOH 

CON'VON  N*RT  (  7  )  tWT  (  7  )  tWR  (  7  )  ?  AP  (  7 , 7  )  .NPf’NT  ,  viyAX  .<MAX  v  ELTA  ♦XTAUt 

1  ZERLAMtXLAMI  2)  »02<7«7),R?(7«7) » I  FLAG  .  P < 2 8 » 1 0 1  )  . T { 149 1 ) . S I G ♦ 

2  P(28»101)»H(28»3»101)  ♦  PT  AU  ♦  Plam  (  2  )  ♦HTAUO)  .HLAM(2*3)  ♦P2{7»7)  ♦ 

3  H2(7.7»3) »CONST(3) »NEQ 
IFLAG=2 
T ( 2 ) »0.0 
T( 3)«DELTA 
M=1 

P'S 

L1=0 

L3  =  3 

DO  1  I=lfN 
DO  1  J=lfl 
Ll'Ll+1 
L3=L3+1 
P{L1»M)»0.0 

1  T(L3)=P(L1  tM,) 

L3=L3+1 
PLAM{2)=0.0 

2  T{L3)=PLAM{2) 

H'S 


00  7  K=l,l 
L1=0 

00  b  I=1»N 
DO  3 
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L3=L3+1 
H{L1  fKjVOsO.O 
3  T(L3)  =H(L1  tlCtM) 

C 

L3=L3+1 

6  HlAM( 2 »< ) = 1 .0 

7  T(L3)=HLAM{24K) 

C 

L  =  0 

DO  8  1=1 fN 
DO  8  J=1»I 
L  =  L+1 

8  R? ( I »J)=R(L»M) 

DO  9  I =1 fN 

DO  9  J-:!  ,N 

9  R2( I ♦J)=R2( J*I  ) 

C 

NEQ  =  2*( (N*(N  +  1 ))/2  +  1  ) 

CALL  INT5(T*NEO»2*0»0»0*0.0»0) 
LMAX=(N*(N+1) )/2 

PRINT52»T(2)*(P(L*M)»H{L»1*K)«L=1 aVAX) 
52  FORMAT(1HOF9.A*5E20.8/{10X5E20.8) } 

C 

DO  51  M1«1»M1MAX 
DO  50  M2=1*NPRNT 
CALL  INTM 
M  =  M+1 

CPREV. APPROX.  R( I tJ) 

L1«0 

DO  10  I-1»N 
DO  10  J*1.I 
L1=L1+1 

10  R2 ( I tJl-RCLl »M) 

DO  11  1*1 tN 

DO  11  J-I.N 

11  R2{  1  tJ)*R2Ui  n 

LI  =0 
L3*3 

DO  12  I*ltN 
DO  12  J*1»I 
L1*L1+1 
L3=L3+1 

12  P(L1  tMl-KLai 
L3=L3+1 

DO  13  K»ltl 
L1=C 

DO  14  I=ltN 
DO  14  J=1.I 
L1=L1+1 
L3=L3+1 

14  H(L1 ♦<»M) *T( L3) 

13  L3=L3+1 

50  CONTINUE 

51  PRINT52»T(2)»(P(L»M).H(L»1»M).L«1 »LVAX  ) 
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pi^TURN 

?ND 

SIBFTC  linear 

SUBROUTINE  LINEAR 
DIMENSION  CHKI ( 3 ) 

DIMENSION  A( 49i3 ) »B (49 ) .EMAT ( 50.50)  .  P I VOT ( 50 ) . I NDEX ( 50 . 2 ) 

l.IPIVOT(50) .FVEC(50.1) 

COMMON  N.RT ( 7 ) . WT ( 7  )  .WR ( 7  )  .  AR ( 7.7  )  .NPRNT .MIMAX .KMAX .DELTA . XTAU  , 

1  Z ERL AM, X LAM ( 2 ),B2(7,7).R2(7.7). I  FLAG .R ( 28 . 101 ).T(1491),SIG» 

2  P(28.101).H(28»3.101 ) . PT AU » PLAM ( 2 ) ,HTAu( 2 ) .HLAM (2.3).P?(7  7)» 

3  H2( 7 .7,3 ) .C0NST( 3 ) •NEO 
CBOUNOARY  CONDITIONS 

MLAST=NPRNT»M1MAX  +  1 
Du  1  K«1 , i 
L  =  0 

00  2  I =1  .N 
DO  2  J  =  1 .  I 
t-  =  L  +  l 

2  H2  (  I  t  J.K)  =H(  L.tC.MLAST  ) 

00  1  I=1.N 

DO  1  J=I.N 

1  H2( I »J.<)=H2( J, I  ,KI 

L  =  C 

DO  3  1=1, N 
DO  3  J  =  ]  ,  I 
L  =  L+1 

3  P2( I»J)=P(L»MLAST) 

DO  4  Ul.N 

DO  4  J=I fN 

4  P2 ( I » J)=P2( J » I  > 

cleast  squares 

DO  5 
L*C 

DO  5  I«1»N 
DO  5  J=1»N 
L  =  L  +  1 

5  A( L.K ) »H? (  I  , J.K ) 

L»0 

DO  6  1=1 »N 
DO  6  J*1»N 
l  =  1 

6  B(L)*82( I .J)  -  P2( I tJ) 

C 

LMAX*N**2 
PR INT60 

60  FORMAT  (,1HQ) 

DO  61  L=1»LMAX 

61  PRINT82. ( A (L»K  )  .K  =  l  .1  )  fB (L ) 

C 

DO  8  1=1.1 
DO  7  J=l.l 
SUM=0.0 
DO  9  L=1.LMAX 
9  SUM=SUM  ♦  A(L.I )*A(L.J) 

7  EMAT( I . J) •SUM 
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SUM=0.0 

DO  20  L=1,LMAX 

10  SUM  =  SUM  +  A(  L»  n#B(L) 

8  FVECn,l)=SUM 

C 

PRrNT60 
00  81  1=1*1 

81  PRINT82» ( EMAT ( I »J) ♦ J«1 ,1 ) *FVEC( I ♦ 1 ) 

82  format { 10X6E20. 8 ) 

C 

FVEC  n  »1 1 »FVEC( 1 *1 ) /EMAT( 1*1) 

C 

00  11  1=1,1 

11  rnw<;T  /  t  \ -c\/rr  »  f  .  *  » 

c 

XLAM<2)=C0NST(1  ) 

PRINT903*XTAU*ZERLAM*XLAM(1) ,XLAM ( 2  ) 

9C3  FORMATdHO/ 

1  IXIIHTHICKNESS  =*  FlO.A  / 

2  1X1 1HALBED0( X )  =*  20HA  *  e>TANH ( 1 0« ( X-C ) )  // 

3  1X3HA  s,  E16.8,  10X3HB  =»  E16.8,  10X3HC  =♦  E16.8  //) 
C 

CNEw  approximation 

r 

M=1 
L  =  0 

00  12  1*1, N 
DO  12  J»1,I 
L  =  L^  1 

SUMsP(l,M) 
do  13  <*1,1 

13  SUM  «SUM  •*  C0NST(<)*H(L»K*M) 

12  R(L,M)*SUM 
L  =  0 

DO  14  I*1*N 
00  14  Jj!  1  ♦  I 
L  =  L  +  1 

14  R2n  ,J)*R(L»M) 

SIG=0.0 

CALL  OUTPUT 

DO  50  M1»1,M1MAX 
DO  18  M2*1*NPRNT 
M  =  M+1 
L  =  0 

DO  15  1=1, N 
DO  15  J=1,I 
L  =  L  +  1 

SUM=P(L,M) 

DO  16  K»l,l 

16  SUM«SUM  ♦  CONST(X.)*H(L,K,M) 

15  R-tLtVlsSUM 
L  =  0 

DO  17  I.1,N 
DO  17  J=l,l 
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l,  =  L  +  l 

17  R2( I »J)»R(L»M) 

18  SIG=5IG  +  delta 
50  call  output 

C 

Rt^TURN 

FND 

tIBFTC  OUTPUT 

SUBROUTINE  OUTPUT 
DIMENSION  X(3) 

COMMON  N»RT(7) ,WT(7)  »WRC')  .AR(7.7>  . NPRNT . M IMAX ♦ KMA X ♦ DELTA , XTAU ♦ 

1  ZERLAM,XLAM( 2 ),B?(7.7).R2(7.7),IFLAG»P(28.1C1)»T{1491)»SIG. 

2  P(28»10n»H(28»3»101)  »  P  T  Au ,  PlAM  (  2  )  .HTAU  I  3  )  ♦HLAM  (  2  *3  )  »P2  (  7  *7  )  ♦ 

3  H2 { 7 3 ) ♦CONST ( 3 ) *NEQ 
DO  1  I=*1»N 

DO  1  J=I»N 
1  R2( I ♦J)=R2( J^I  ) 

Y=XTAU*SIG 
X ( 1 JsZERLAM 
X ( 2 ) =XLAM{ 1 ) 

X( 3)*XLAM( 2) 

CALL  ALBEDO! Y,X,Z) 

PRINTIOO,  SIG.Y»Z 

100  FORMATdHO  7HSIGMA  =tF6»2*  4X5HTAU  =♦  F6.2*  4X8HALBEDO  =^F6.2/) 
DO  2  J=lfN 

2  PRINT101.J.(R2{ dJ)  ♦I'lfNJ 

101  FORMATdlO^  7F10,6) 

RETURN 

END 

SIBFTC  ALBEDO 

SUBROUTINE  ALBEDO ( Y  ♦  X  ,  Z ) 

DIMENSION  X(35 

COMMON  NfRT(7) .WT(7i  .WR(7)  fARITfT)  ♦ NPPNT ♦M 1 MAX ♦ KMA X ♦ DEL TA ♦ X TAU ♦ 

1  ZERLAM,XLAM( 2) ♦B2(7f7).R2(7f7),IFLAG^R{28tlCl)^T(149l)>SIG> 

2  P(28.101)fH(28>3fl01) tPTAUtPLA'^  ( 2 ) ♦HTAU( 3 ) ♦HLAM( 2  ^3 )  ♦P2 ( 7,7)  ♦ 

3  H2{7,7,3)  ♦CONSTO)  ♦NEQ 
ARG=10.0*( Y-X ( 3  )  ) 

Z  =  X( 1  )  +  X(2)*TANh(ARG) 

RETURN 

END 
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11^  .program  for  the  estimation  of  t,  the  thickness  of  the  medium 
SIBFTC  RTINV 

COMMON  N»kT( 7)  ,WT(7)  »WR(7)  >AR(7«7 J  . NPRNT , M 1 MAX t <M A X  ♦  DEL TA  ♦  X T AU 

1  Z ERL AM, XL AM ( 2 ),B2(7»7!»R2(7»7», IFLAG»R( 28 ♦ 1 01 ) ♦ T (149 1 ) ♦ S I G » 

2  P(28,101!»H(28,3»101! , P T AU , PLAM ( 2 » ,HTAU(3) .HLAM(2»3! ♦P2(7,7) 

3  H2(7,7,3! ,C0NST(3).NF0 

C 

C  PHASE  I 

C 

1  READ1000,N 
PRINT899 
PR  INT900  tN 

RFADlOOl , ( RT ( n  ,  1  =  1  ,N) 

PRINT901 ♦ (RT( n  ,  1  =  1  ,N) 

READlCOl, (WT( I ) ,I=1,N) 

PRINT9ni,(WT'  I J  »I  =  1»N) 

DO  2  1=1»N 

WR( I !=WT( I )/RT( I ) 

DO  2  J=1,N 

2  AR(I,J)=  l.O/RTd)  +  1,0/RT(J) 

C 

899  F0RMAT{1H146XA-HRADIATIVE  TRANSFER  -  INVERSE  PROBLEM  /  ) 

1000  F0RMAT(6I12) 

900  FORMAT(6I20) 

1001  format (6E12. 8 ) 

901  FORMAT (6E20. 8 ) 

R EAD 1 000, NPRNT, Ml MAX, KMAX 
PRINT900,NPRNT,M1MAX»KMAX 
READlOOl, DELTA 
PRINT901, DELTA 

READU'Ol ,XTAU,2ERLAM»XLAM{ 1 ) ,XLAM(2) 

PRINT902 

PRINT903.XTAU»2ERLAM»XLAM{ 1 ) ,XLAM(2) 

902  FORMAT! 1H123HPHASE  I  -  TRUE  SOLUTION  /) 

903  FORMAT (IHO/ 

1  IXIIHTHICKNESS  »»  FlO.4  / 

2  1X11HALB£D0(X)  =,  20HA  +  B*TANH ( 10* ( X-C J )  // 

3  1X3HA  »,  E16.8»  10X3HB  E16,8,  10X3HC  “*  E16.8  //) 

CALL  NONLIN 

DO  3  1=1, N 
00  3  J«1,N 

3  B2 ( I tJlsRZI I tJ) 

C 

C 

C  PHASE  II 

C 

4  REA01001,XTAU,ZERLAM,XLAM( 1 ) ,XLAM{2» 

K»C 

PRINT904,K 

PRINT903,XTAU,ZERLAM,XLAM( 1 ) ,XLAM(2) 

C 

CALL  NONLIN 
C 

904  FORMATdHl  13HAPPR0X  IMAT  I  ON ,  13/  ) 

C 

C  OUASILINEARIZATION  ITERATIONS 


V-'  KJ 
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DC  5  Kir  1 ,KMAX 
ORINT<>04»K1 
CALL  PANDH 
CALL  LINEAR 
5  CONTINUE 


C 

READlOOO* IGO 
GO  TO  { 1 *4) ♦ IGO 
END 

SlBETr  HAUX 

subroutine  DAUX 

DIVFNSION  V2(7»7),X(35  ♦F(7)»G(7) 

CO^^MON  N»RT(7)  ♦WT(7)  tWRl?)  fARtT.^)  ,  NP^NT  ..MIMAX  tKMAX  tDELTA  ♦  XTAU  ♦ 

1  ZERLA,-1»XLAM(  2  )  ♦82(7»7)»R2(7»?).  IF  LAG. P{28»101)  ♦1(1491  )»5IG» 

2  9(28.101  ) 28 »3» 101) . PTAU . PLAM ( 2 ) .HTAU( 3 )  .HLAM ( 2 ♦ 3 )  ♦  P2 ( 7 . 7  )  . 

3  H2(7.7.3) »CONST(3 ) ♦NEQ 
GO  TO  (1.2) » IFLAG 

C 

CNONLINEAR 

r 

-w 

.  L  =  3 

DO  4  1=1 .N 
DO  4  Jal.I 
L  =  L+1 

4  V2(I»J)=T(L) 

DO  5  I=1.N 
DO  5  J=1 .N 

5  V2( I »J)*V2( J«I  ) 

L  =  L  +  1 

VT^U«T( .f 
SIG»T(2) 

YsVTAU-sSIG 
X(  1  I-IEPLAM 
X(2)«XLAM( 1 ) 

X(7).XLAM(?) 
call  albedo* Y.X.Z) 

ZLAMDAaZ 

C 

DO  6  I»1 .N 
F(  I  )--0.0 
DO  7  K«1.N 

7  F(  I  )«F( I)  ♦  WR(K)*V2( I .K) 

6  F( i )a0.5*F( I  )  4  1*0 
C 

Du  8  I«1*N 
DO  e  J«l»l 
L  =  L  +  1 

DR«-AR!  I  .J)»V2(  I  tJ)  ♦  ZLAMOA*.'M  n  *F  ( 

8  T(L)*OR*\/TAU 
DO  9  I-l.l 
L-L  +  l 
T(L)»0,0 
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o 


RETURN 
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CLINEAR 

C 

2  SIG=T*2) 

Y=XTAU*SIG 
X(  1  )=ZERlAM 
X ( 2  )  =XLAM( 1  ) 

X ( 3 ) *XLAM ( 2 ^ 

CALL  ALBEDO(V.XfZ) 

ZLAMDArZ 

C 

DO  16  I«1»N 
F(  I  )»0.0 
DO  17  K=1,N 

17  F(I)=F(I)  •>  WR(IC)*R2(  I  »IC) 

15  F( I )*0.5*F( I )  ♦  1.0 
C 

CP'S 

C 

L  =  3 

DO  14  1=1, N 
DO  14  J*1,I 
L  =  L  +  1 

14  V2(I,J)=T(L) 

DO  15  J=1,N 
DO  15  J= I ♦ N 

15  V2( I ,J)»V2( J»  I  ) 

L  =  L  +  1 
VTAU«T{L) 

C 

DQ  10  1  =  1, N 
G{  I  )»0.0 
DO  10  K=1,N 

10  G{I)=G{I)  +  {V2{  I  ♦K)“R2n  ♦<)  )*WR(IC) 

ARG=1C.0*( Y-XLAM{2) J 

PARTL  =  10.0*SrG*XLAM{ 1 ) ♦( 1 . 0- ( TANH ( ARG )  )  **2  ) 
M=3+NEO 
DO  12  1=1, N 
DO  12  J»1,I 
FIJ=F{ I )*F( Jj 

CAPF  =  -mR(  I  ,  Jl*R2  {  I  »  J)  ZLAMOA*FIJ 
Tla-XT; J*AR{ I ,J)*( V2{ I f J)“R2{ I t J)  ) 
T2=^XTAU*0.5*ZLAMDA*{F{  I  )  *G  {  J  )  ♦F  !  J  )  *rM'  I  )  ) 
T3=XTAU*CAPF 

T4=  (  V!'AU-XTAU)*(CA'5F  +  XT  AU*F  IJ»P  ARTL  ) 

M=M+1 

12  T(M)=T1+T2+T3+T4 

DO  19  I«l,l 
M  =  M+1 

19  T(M)*0,0 

C 

CH'S 

C 


DO  100  K.l  ,1 
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C 

00  24  I»1.N  c 

00  24  J»1,I 
L  =  L  +  1 

24^  V2(I»J)*T(L)  C 

DO  25  I«l.N 
00  25  J»I .N 

25  V2(  I  .J)»V2(  J*n  C 

L‘L+1 
VTAU'T(L) 

C 

DO  20  I»1»N 
G( I 1»0.0 
DO  20  J*1*N 

20  G(:)»G(I)  -f  V2n»J)*WR(J) 

C 

00  22  I»1»N 
00  22  J»1»I 
FIJ=F(I )#F( J) 

CAPF»“AR( I , J) *R2( I ♦ J)  ♦  ZLAMDA*FIJ 
T1»“XTAU*AR( I ♦J)*V2 ( I » J) 

T2»XTAU*0.5*ZLAM0A*<F( I ) *G ( J J +F ( J ) *0 ( I ) ) 

T3»0.0  C 

T4*VTAU»(CAPF  XTAU*F  IJ*PARTL) 

V  =  M+1 

22  T (M) •T1+T2+T3+T4  S 

C 

00  29  I*l»l 
M«M+1 

29  T(M)»0*0 

100  CONTINUE 

RETURN 
END 

SIBFTC  NONLIN 

SUBROUTINE  NONLIN  Cl 

COMMON  N»RT(7) ,WT(7) »WR(7) fAR(7»7) , NPRNT , Ml MAX . KM AX , DEL TA , XT AU » 

1  ZERLAM»XLAM( 2 ) .02 17»7 J ,R2 (7»7) , I FL AG  .9 ( 28 » 1 0 1 ) » T ( 149 1 ) . S I G . 

2  P(28»101)»H(28»3.i01) »PTAU.PLAM( 2 ) .HTAU( 3) .HLAMI 2*3'  .P2 ( 7,7  )  » 

3  H2(7,7,3) »C0NST(3 ) ,NEQ 

C  NONLINEAR  D.E.  FOR  TRUE  SOLUTION  OR  FOR  INITIAL  APPROX, 

C 

IFLAG'l 

T{2)»0.0 

T(3)«0ELTA 

Mel 

L1»0 
L3  =  3 

00  1  1*1. N 
00  1  J=1,I 
ll=Ll+l 
L3«=L3+1 
R2 ( I .J)=0.0 
R(L1,M)=R2( I .J) 

1  T(L3)=R2n,J)  Cl 
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L3=L3+1 

2  T(L3)=XTAU 
C 

NE0=(N»(N+1 ) )/2  +  1 

call  !NTS(TtNEOt2tOtOtOfOtO»0) 

C 

5IG  =  T (2  5 
CALL  Output 
c 

DO  5  M1*1,M1MAX 
DO  4  M2=1»NPRNT 
CALL  INTM 
M  =  M+1 
L1=0 
L3  =  3 

DO  3  1=1. N 
DO  3  J=1.I 
Ll=Ll-»  1 
L3  =  L3-^1 
R2(I.J)»T(L3) 

3  R(L1 .M)*R2( I .J) 

4  SIG=T(2) 

5  CALL  OUTPUT 
C 

RETURN 

END 

SI0FTC  LINEAR 

SUBROUTINE  LINEAR 
DIMENSION  CHKI ( 3) 

DIMENSION  A{49.3)  .8(49) .EMATISO.SO) .  P I VOT ( 50 ) . INDEX { 50 . 2 ) 

l.IPIVOT(50) .FVEC(50.1) 

COMMON  N.RTt  7) .WT{7) .WR{7) .AR ( 7 . 7 ) . NPRNT .MIMAX . KMAX .DEL TA . XTAU . 

1  ZERLAM.XLAM( 2) .82(7.7) . R2 ( 7 , 7 ) . I  FLAG .R ( 28 . 101 ) .T ( 1491 ) .SIG. 

2  P(28.101 ) »H(28.3. 101 ) .PTA'J.PLAM ( 2 ) .HTAU( 3) .HLAM(2»3) *P2( 7.7) * 

3  H2(7.7.3).CONST(3).NEO 
CBOUNDARY  CONDITIONS 

MLAST=NPRNT*M1MAX  ♦  1 

DO  1  K*l.l 

L«0 

DO  2  I=1.N 
DO  2  J*1.I 
L  =  L  +  1 

2  H2(I.J.K)«H(L.)C.MLAST) 

DO  1  1=1. N 

DO  1  J'l.N 

1  rf2  (  I  .J.IO«H2(  J.  I  .K) 

L  =  0 

DO  3  I«1.N 
DO  3  J«1.I 
L'L-*-! 

3  P2(I.J)=P(L»MLAST) 

DO  4  I*1.N 

DO  4  J=I.N 

4  P2{ I »J)*P2(J.I  ) 

cleast  squares 
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DO  5  K=l.l 
L  =  0 

DO  5  I=1.N 
DO  5  J=1»N 
L  =  L  +  1 

5  A{L»<)=H2( I tJtK) 

L  =  0 

DO  6  1=1. N 
DO  6  J=1.N 
L  =  L+1 

6  B(L)=B2( I . J)  -  P2( I .J) 

LMAX=N**2 

PRINT60 

60  FORMAT (IHO) 

DO  61  L=1.LMAX 

61  PRINT82. (A(L.K) .K«l »1) .B(L) 

DC  8  1=1.1 
DO  7  J=l.l 
SUM=0,0 
DO  9  L=1.LMAX 

9  SUM=SUM  +  A(L.I )*A(L.J) 

7  EMAT(  I..J)  =SUM 
SUM=0.0 

DO  10  L=1.LMAX 

10  SUM=SUM  +  A(L.I)*8(L) 

8  FVEC( I .1 ) «SUM 

c 

PRINT60 
DO  81  1=1.1 

81  PRINT82.(EMAT{I.J).J»1.1).FVEC(I.1) 

82  FORMAT(10X6£20,8) 

C 

FVEC ( 1 ♦ 1 ) -FVEC ( 1 ♦ 1 ) /EMAT (1.1) 

C 

DO  11  1=1.1 

11  CONST( I )«FVEC( I  .1 ) 

C 

XTAU  «CONST(l) 

PRINT903»XTAU»ZERLAM.XLAM( 1 ) .XLAM(2  ) 

903  FORMATdHO/ 

1  IXIIHTHICXNESS  =.  E16.8  / 

2  IXllHALBEOOIX)  «.  20HA  ♦  B*TANH ( 10*( X-C) )  // 

3  1X3HA  =♦  E16.8.  10X3HB  =.  E16.8.  10X3HC  =♦  E16.8  //) 
C 

CNEW  APPROXIMATION 
C 

M»  1 
L  =  0 

DO  12  1=1. N 
DO  12  J*1.I 
L  =  L*1 

SUM*P(L  »M) 

DO  13  <=1.1 


%j  KJ 
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13  SUM  =SUM  +  C0N5T{K)*H{L»K.M) 

12  RCL»M)=SUM 

L  =  0 

00  14  Ul.N 
DO  14  J2i,l 
L  =  L+1 

14  R2 { I » J)*R ( L»M» 

SI6=0.0 

CALL  OUTPUT 
C 

DO  50  M1*1,M1MAX 
DO  1.8  M2=1«NPRNT 
M  =  M+1 
L  =  0 

DO  15  I«1»N 
DO  15  J=1.I 
L  =  L  +  1 

SUM=P (LiM) 

DO  16  <=1»1 

16  SUM=SUM  +  CONST(K)#H(L»< fM) 

15  R{L.M)=SUM 
L  =  0 

DO  17  I»1,N 
DO  17  J*1,I 
L  =  L  +  1 

17  R2 ( I » J)«R ( L»M) 

18  SIG=SIG  +  DELTA 
50  CALL  OUTPUT 

C 

RETURN 

END 

SIBFTC  PANOH  list 

SUBROUTINE  PANDH 

COMMON  N»RT ( 7  )  »WT ( 7 ) »WR ( 7 ) »AR ( 7»7 ) fNPPNT  ^MIMAX  »KMAX  »DELTA ,X  TAU 

1  Z£RLAM,XLAM(2) »B2(7.7) . R2 ( 7 ♦ 7 ) , I  FLAG ,R ( 28 » lOl ) »T{ 1491 ) ,SIG, 

2  P(28»101)»H(28»3»101) »PTAU»PLAM { 2 ) »HTAU{3) »HLAM { 2 »3 ) ♦P2 ( 7 »7  ) 

3  H2(7,7»3)»CONST(3)»Nen 
IFLA^=2 
T(2)«0.0 
T( 3)»DELTA 
M»1 

P‘S 

Ll-0 
L3  =  3 

DO  1  1=1, N 
DO  1  J=1,I 
L1*L1+1 
L3»L3+1 
P(L1 ,M)*0.0 

1  T(L3) =P(L1 fM) 

L3=L3+1 
PTAU*0.0 

2  T(L3)=PTAU 
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C  M*S 
C 

DO  7  k;*i.i 
L1=0 

DO  3  I=liN 
DO  3  J=1»I 
L1=L1+1 
L3  =  L3-*-l 
H(Ll»K.M'j=0.0 
3  T(L3)=H{L1 .K»M) 

C 

L3=L3+1 

6  HTAU(K)=1.0 

7  T(L3)=HTAU(K) 

C 

L=0 

DO  8  1=1 .N 
DO  8  J=1 ♦ I 
L  =  L  +  1 

8  R2( I »J)=R(L»M> 

DO  9  I=1»N 

DO  9  J=I»N 

9  R2 ( I  ♦  J  )  =  R2 ( J* I ) 

C 

NEG=2*(  (N*(N+1  )  )/2  1  ) 

CALL  INTS(T»NEQ*2*0*0»0»0»0»0) 
LMAX=(N*(N+1 ) )/2 
C 

DO  51  M1=1,M1MAX 
DO  50  M2*1»NPRNT 
CALL  INTM 
M  =  M+1 

CPREV. APPROX.  R( I »J) 

LI  *0 

DO  10  1*1  ,N 
DO  10  J»1,I 
L1=L1+1 

10  R2n  *J)»R{L1»M) 

DO  11  I»1 .N 

DO  11  J=I fN 

11  R2 (I »J)=R2{ Jt I ) 

L1=0 

L3*3 

DC  12  1*1. N 
DO  12  J*1 . I 
L1*L1+1 
L3=L3+1 

12  P(L1.M)=T{L3) 

L3=L3+1 

DO  13  K=l.l 
L1=C 

DO  14  1=1. N 
DO  14  J=1»I 
L1=L1+1 
L3*L3+1 
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lU  H ( L 1 »K »M) =T ( L3 ) 
n  L'’=L''  +  1 

50  rONTiNUF 

51  rONTINU'^ 
return 
END 

sierrc  output 

subroutine  Output 
dimension  X{3> 

COMMON  N»RT ( 7 ) »WT  I  7 ) .WR  (  7 ) .AR ( 7»7 ) »NPRNT ,M1MAX  »<MAX tOELTA  tXTAU  ♦ 

1  ZERLAM,XLAM,( 2 ).B2(7.7).R2(7»7)»IFLAG,R(28.101)»T(1491).SIG» 

2  P(28»101)»H(28»3.101) . P T AU . PLAM ( 2 ) tHTAUl 3) tHLAMI 2  *3 ) ♦P2 ( 7.7  )  ♦ 

^  H2(7.7.3) *CONST(3) »NEO 

DO  1  I=1.N 
DO  1  J=I.N 
1  R2( I »J)=R2( J» 1 ) 

Y=XTAU*SIG 
X( 1 ) =ZERLAM 
X ( 2 ) =XLAM(  1  ) 

X( 3)=XLAM( 2) 

CALL  ALBEDO(Y.X.Z) 

PR  I  NT  100.  SIG.Y.Z 

100  FORMATdHO  7HSIGMA  =.F6.2»  4X5HTAU  =.  F6.2.  AX8HALBEDO  =.F6.2/) 
DO  2  J=1.N 

2  PRINT101.J.(R2( I .J)  .I  =  1.N) 

101  FORMATdlO.  7F10.6) 

RETURN 

END 

5IBFTC  ALBEDO 

SUBROUTINE  ALBEDO ( Y  .  X  .  Z  ) 

DIMENSION  X(3) 

COMMON  N.RT( 7) .WT( 7)  .WR(7)  .AR(7.7)  .  NPRNT .Ml  MAX . <MAX . DEL TA . XTAU  . 

1  ZERLAM.XLAM(2)  .82 ( 7,7)  .R2 (7.7) . I  FLAG .R ( 20 . 1 0 1 ) . T ( 149 1 ) ♦ S I G  . 

2  Pv23»101)»H(28.3»101) .PTAU .PLAM ( 2 ) .HTAU ( 3 ) »HLAM(2»3) .P2(7.7) , 

3  H2(7.7.3) »C0NST(3) ♦NEO 
ARG*10.0# ( Y-X ( 3  )  ) 

Z  =  X( 1 )  ♦  X(  2)»TANH{ ARG) 

RETURN 

END 
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ni; _ PROGRAM  FOR  THE  ESTIKATIQN  OF  X,,  X^,  c,  THE  TWO  ALBEDOS  AND  THE  THICKNESS 

^  38*^  OF  THE  LOWER  LAYER 

SIBFir  RTINV 

COMMON  N , RT  <  7  )  » WT  <  7  )  »WR  <  7  )  » AR I ?*  7  )  ,NORNT  *M1MAX  »KMAX  *DELTA  ,  XTAU  ♦ 

]  XLAM{3).  B2 <7,7)  ,R2(7»7) , lFLAG»R(28»ini ) ,T( 1491 )  .SIG* 

?  P(28»l01)»H<28»8.101)»PLAMi3) *HLAy(3»3) iP2(7»7) . 

3  H2  (  7 , 7, 3- )  ’CONST  <  3  )  »N£Q 

phase  I 

1  READIOOO’N 
PRINT899 
PRINT900.N 

RFADlOOl ♦ ( RT ( I  )  »I  =  1  »N) 

PRINT901, (RT( I)  »!  =  1  ’N) 

RFADlCOl » ( WT (  I  )  ,  1  =  1 »N) 

P'RINT901»  (WT(  I  )  »I  =  1  »N) 

DO  2  I=1»N 

WR( I ) =WT( I )/RT( I ) 

DO  2  J=1»N 

2  AR(I»J)=  1.0/RT(I)  1.0/RT(J) 

899  format ( IH]  46X36HRA0 I  AT  I VE  TRANSFER  -  INVERSE  PROBLEM  /  ) 

1000  F0RMAT(6I12) 

900  FORVAT(6I20) 

1001  format (6F12. 8 ) 

901  FORMAT(6F20.8) 

REAOl  OOOtNPRNT  »M1MAX  »<MAX 
PR lNT90CtNPRNT  tMlMAXtKMAX 
REA.OlOOl  ’DELTA 
PRINT901»DELTA 

READlLOl ’XTAU* ( XLAM( I )  » I^l t3) 

PR  INT902 

PRINT903»XTAU»(XlAM( I ) »I«1»3) 

902  FORMAT ( 1H123HPHASE  I  -  TRUE  SOLUTION  /) 

9C3  FORMAK  IHO/ 

1  1X1 IHTHICKNESS  «»  FlO.4  / 

2  I'Xl  1HAL8E00(  X )  «»  20HA  +  0*TANH(10*(X-C)  )  // 

3  1X3HA  »♦  E16.6’  10X3HB  =« »  E16.8»  I^x^hC  =>  E16.8  //) 

call  nonlin 

DO  3  I=1’N 
DO  3  J=liN 

3  B2n  t J)«R2<  I  fJ) 


PHASE  II 

4  READ1001.XTAU»(XLAM(  n  »I«1*3) 

K=r 

PRINT904»K 

PRINT903»XTAU» ( XLAM( 1) »I»1»3) 

C 

CALL  NONLIN 
C 

904  FORMAT<lHi  1 3HaPPROX I  MAT  I  ON »  13/  ) 


C 

C 


QUASILINEARIZATION  ITERATIONS 


n  o  r> 


C 
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DO  5 

PR  I NT904»K1 
CALL  PANDH 
CALL  LINEAR 
5  CONTINUE 


READIUOO, IGO 
GO  TO  ( 1 »A ) » IGO 
END 

SIBFTC  DAUX  LIST 

SUBROUTINE  DAUX 

DIMENSION  V2I 7,7) »X ( 3 ) tF ( 7 ) fGI 7 ) 

1  »VLAM(3) 

COMMON  N,RT(  7)  ,V»T  I  7  )  fWRI  7)  ,AR  (  7, 7  )  .  NPRNT  ,  M IMAX  » KMAX  »OEL  TA  ,  XT  AU  » 

1  XLAM(3),  B2<7,7) ,R2(7,7) , IFLAG,R(28»101 ) *T(1491) ,SIG* 

2  P(28,101)»H(28,3,101),PLAM(3) ,HLAM( 3 ,3)  ,P2{ 7,7)  , 

3  H2(7,7,3) »CONST(3)»NEQ 
GO  TO  (1,2)  ♦  IFlAG 

C 

CNONLINEAR 

c 

1  L  =  3 

DO  A  I=1,N 
DO  A  J=1,I 
L  =  L  +  1 

A  V2(I»J>=T(L) 

DO  5  I=1,N 
DO  5  J=I,N 

5  V2 ( I t J)=V2 ( I  ) 

DO  51  I»l,3 

L  =  L  +  1 

51  VLAM{I)=T(L) 

SIG=T(2) 

Y=XTAU*SIG 
DO  52  1=1,3 

52  X(  I  ):fVLAM{  I  ) 

CALL  ALBEDO! Y,X,Z) 

zlamda=z 

C 

DO  6  I =1  ,N 
F ( I ) =0.0 
DO  7  K=1,N 

7  F( I  )*F( I)  +  WR(<)*V2( I  ♦<) 

6  F( I  )»C.5*F(  I  )  +  1.0 
C 

DO  8  I=1,N 
DO  8  J=1,I 
L  =  L  +  1 

0R=-AR I  I , J )*V2( I t J)  +  ZLAMOA»F( I )*E(J) 

8  T(L)=DR 

DO  9  1=1,3 
L  =  L  +  1 


. . 
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9  T(L)=C.O 
RETURN 
C 

CLINEAR 

C 

2  EfG=T(2) 
y=xTAU*SIG 
DO  21  1=1 .A 
21  X ( I ) =XLAM( I ) 

call  ALRED0( Y.X.Z) 

2LAyDA=Z 

C 

DO  16  1=1. N 
F(  1)  =0.0 
DO  17  <=1.N 

17  F(  I  )=F(  I  )  >  WR(<)»R2{ I »<) 

16  F( I )=0.5*F( 1 )  +  1.0 

C 

CP'S 

c 

L  =  '» 

DO  14  I=l.N 
DC  14  J=1 ,  I 
L  =  L>1 

14  V2(I.J)»T(L' 

DO  15  I=1.N 
DO  15  J=I.N 

15  V2( I ♦J)»V2( I ) 

DO  18  1*1.3 
L*L+1 

18  VLAM(I)»T(L) 

C 

DO  10  1=1. N 
G(  I  )*0,n 
DO  10  <=1.N 

10  G(I)=G(I)  +  (V2( I»K)-R2n .<) )*WR(K) 

ARGsl0.0*(Y-XLAM(3)  ) 

TARG=TANH( ARG) 

XTANX=-10.0*XLAM( 2)*( 1 .0“TARG**2 ) 
M=3+NE0 
DO  12  1  =  1  .N 
DO  12  J=1.I 
FIJ=F(I)*F( J1 

CAPF  =  -AR{  I  .  J)  »R2(  I  .  J)  ZLAMDA*FIJ 
T1 =C4PF 

T2=-AR ( I . J )* ( V2 ( I . J)-R2(  I  ♦  J  )  ) 

1  +  0.5*ZLAMDA*(F(  I  )*G(  J)  F(J)*G(I)) 

T3=(VLAM( 1 )-XLAV( 1 )  )*FT J 
T4=  (  VLA.V(  2  ) -XLAM(  2  )  )*TARG*FIJ 
T5»(VLAM(3)-XLAM(3) )*XTANX*FIJ 
M  =  M+1 

12  T (M ) =T1+T2+T3+T4+T5 
DO  19  1=1.3 
M  =  M+1 
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19  T ( M ) =  0 . 0 
C 

CH*  s 

c 

DO  lOC  <a  1 3 

c 

00  ?4  1*1 tN 
DO  24  J=1,I 
L  =  L-t-l 

24  V2(I»J)=T(L) 

DO  25  1*1 tN 
DO  25  J= 1  .N 

25  V2 ( I ♦ J)=V2 ( J»  I) 

DO  26  I=l»3 

L  =  L+1 

26  VLAM{I)*T(L) 

C 

DO  20  I=1,N 
G( I )=0.0 
DC  20  J=1,N 

20  G(I)»G(I)  +  V2( I vJ)*WR( J) 

C 

00  22  1*1  . N 
DO  22  J«1»I 
FIJ=F( I )*F( J) 

T1»0.U 

T2*-AR(  I  |J)*V2(  I  tJ)  0.5*ZLAMr)A*(F(  n*G(  J)  4  F(J)#G{I)) 
T3*VLAM(1)#F!J 
T4«VLAM(2)*TARG*FIJ 
T5*VLAM(3)*XTANX*FIJ 
M«M+1 

22  T(M)«Tl4T24T3'*-T44T5 

C 

00  29  I«l,3 
M  =  M+1 

29  T(M)»0,0 

100  CONTINUE 
RETURN 
FNO 

SIBFTC  NONLIN 

SUBROUTINE  NONLIN 

COMMON  N,RT( 7) ,WT(7) fWR(7) »AR(7»7) ♦NPRNT,M1MAX»<MAX,DELTA»XTAU» 

1  XLAM{3)»  B2(7»7) ♦R2(7»7) ♦ IFLAG»R(28. 101 ) tTI 1491 ) tSIG* 

2  P(28»101)»H(28f3il01)»PLAM(3)»HLAM(3f3) »P2(7»7)  i 

3  H2(7,7f 3) ♦C0NST(3)iNE0 

C  NONLINEAR  O.E.  FOR  TRUE  SOLUTION  OR  FOR  INITIAL  APPROX. 

C 

IFLAG»1 

T( 2)«0,0 

T(3)*DELTA 

M*1 

L1=0 

L3  =  3 

DO  1  I«1»N 
DO  1  J«1»I 


tNi  cn 
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LI =L1+1 
L3=L3+l 
K  ?  (  I  » J  )  =  0 . 0 
ftCLl »M)=R2( I fJ) 

1  T{L3)=R2(I»J) 

00  2  I=I»3 
L3=L3+1 

2  T(L3)=XLAM( I  ) 

C 

NEO=(N*(N+l )  )/2  +  3 

CALL  INTS(T»NEQ*2*0»0t0»0*0»0) 

C 

SIG*T(2) 

CALL  OUTPUT 
C 

no  5  M1=1,M1MAX 
DO  4  M2=1»NPRNT 
CALL  INTM 
M  =  M+1 
LI  -'0 
L3  =  3 

DO  3  I=1»N 
DO  .3  J  =  1»I 
L1=L1+1 
L3  =  L3  +  ?. 

R2  i  I  »J)sT(L3  ) 

3  R (LI  tMlrRZf I »J) 

4  SIG*T.2) 

■3  call  output 
C 

RETURN 

END 

iI3PT<'  PANDH 

SUBROUTINE  '^ANL  ! 

COMMON  N»RT  (  7)  »Wr  (  7  )  »WR  (  7  )  »  AR  (  7»7)  tNPRNT  ,M1  MAX.  ♦  KMAX  » DELTA  »  XT  AU  ♦ 
1  XLAM(3)»  t2(7»7) ,R2(7»7)»IFLAG»R(2e.l0l).T(149l)»SIG» 

P(23»101)»H(28»3»1C1)»PLAM(3) ♦HLAvi 3,3) ♦P2(7»7) * 

H2(7»7»3) *C0NST(?I »NE0 
IFLAG»2 
T(2)*0.0 
T( 3)«DElTA 
M*  1 

C  P'S 

c 

LI  *0 
L3a3 

DO  1  I»1>N 
DO  1  J=lfl 
L1»L1+1 
L3=L3+1 
P(L1 »M)sO.O 
1  T(L3)=P<L1»M) 

DO  2  i*:»3 
L3*L3+1 
E'.AMI  I  )»0*0 


r»  r> 
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2  T(L3)=PLAM(  {  ) 

H*S 

DO  7  K=l,3 

L1=0 

DO  3  I=1.N 
DO  3  J  =  1 .  I 
L1=L1+1 
L3=L3+1 
H(L1 »K.M) =0.0 

3  T(L3)=H(L1  . 

DO  7  1=1,3 
L3=L3+1 
HLAM( I ,K) =0.0 
IF( 1-K)7t6,7 
6  HLAM(  I  ,K) =1 .0 
'7  T(L3)=HLAM( 1 ,K1 

L  =  0 

DO  8  1  =  1. N 
DO  8  J  =  1  . 1 
L  =  L+1 

8  R2( I »J)=R(L»M) 

DO  9  1  =  1  ,N 
DO  9  J=1,N 

9  R2 ( I ♦ J)=R2 ( J»  1  ) 

NE0  =  4»{  (N*«N  +  1)  )/2  3) 

CALL  INTSl T,NEO,2,0»0,0»0,0»0) 
LMAX=(N#(N+1)  )/2 

RFINT52»T(2)»tP(L»M)»H(L,lfM),L=l »LMAX ) 
52  FORMAK 1H0F9.4,5E20.8/ ( 10X5E20.8) ) 

DO  51  M1=1,M1MAX 
DO  50  M2=1»NPRNT 

call  INTM 

M  =  M+1 

CPRcV. APPROX.  R : I , J) 

L1=0 

DO  10  1  =  1, N 
DO  10  J=1,I 
L1=L1+1 

10  R2( I »J)=R(L1.M) 

DO  11  I=1,N 

DO  11  J*I,N 

11  R2  :  I  ♦  ')=R2(  J»  I  ) 

L1=0 

L?  =  3 

DO  12  1  =  1, N 
DO  12  J=1,I 
LI ^Ll+1 
L3=L3+1 

12  P(L1 .M)=T( L3) 
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L3=L3+3 
DO  13  K=l»3 
L1=0 

DO  14  I»1,N 
DO  14  J=1,I 
Li =L1+1 
L3=L3+1 

14  HfLl L3) 

13  L3=L3-3 

50  CONTINUE 

51  PRINT52tT(2)*(P;L»M),H(L*l*M).L=l  .L^^AX  ) 

RETURN 

END 

5IBFTC  linear 

subroutine  linear 

DIMENSION  CHKI(3) 

DIMENSION  Al 49,3) *8 (49 ) ,EMA1  ( 50»50)  ♦  P I VOT ( 50 )  » I NDEX { 50 ♦ 2  ) 

1 ♦ IPIVGT(50) ♦FVEC(50»1) 

COMMON  N»RT ( 7  )  ,WT ( 7  )  »WR( 7 ) tAR ( 7 .7  )  ,  NPPNT , M 1  MAX  » KMA X ♦ DEL TA , XT AU  ♦ 

1  XLAMO),  B2(7»7)  ♦R2(7.7)  »  IFLAG»R{28»1Q1  )  *7  (  1491  )  .SIG. 

2  P(28»10n»H(28.3»101)  fPLAMI  3  )  .HLAMO  ,3)  ♦82(7.7), 

3  H2 ( 7 .7.3 ) .CONST ( 3 ) .NEC 
CeOUNDARY  CONDITIONS 

MLAST=NPRNT*M1MAX  +  1 
DO  1  K=1.3 
L  =  0 

DO  2  I»1.N 
DO  2  J»1.I 
L  =  L'*-1 

2  H2( I »JtK)sH(L.K.MLAST) 

DO  1  1*1, N 

DO  1  J»I ,N 

1  H2( I ♦J.<)*H2( J,I ,<) 

L  =  0 

DO  3  i»l.N 
DO  -3  J=1 .  I 
L  =  L+1 

3  P2( I »J)*P(L.MLAST) 

DO  4  '=1.N 

DO  4  J=I,N 

4  P2( I ,J)»P2( J» I ) 

CLEAST  SQUARES 

DO  5  <*1.3 
L  =  0 

DO  5  I*1,N 
DO  5  J*1,N 
L»L+1 

5  A(L.K)=H2(  I  . J,tC) 

l=c 

DO  6  U*1,N 
DO  6  J*1,N 
L*L-*-l 

6  B(L)*B2( I tJ)  -  P2( I . J) 

C 

LMAX»N**2 


45 


PRINT60 

60  format (IHO) 

DO  61  L=1.LMAX 

61  PRINT82.<A(L.K)  fK:=1.3)  *0(1) 

C 

DO  8  I=l»3 
DO  7  J=l,3 
SUM=0.0 
DO  9  L=1»LMAX 

9  SUM=SUM  +  A(L. I )*A(L»J) 

7  EMAT( I .J) =SUM 
SUM^O.O 

DO  10  L=1.LMAX 

10  SUMsSUM  +  A(L.I 

8  FVECI I il)«SUM 
C 

PRINT60 
DO  81  1  =  1 »  3 

81  PR  I  NT 82 » ( EMAT ( I « J ) » J= 1 « 3 ) « F VEC M ♦  1  ) 

82  FORMAT ( 10X6E20. 8 ) 

C 

call  MATINV(EMAT«3»FVECf 1 tDETERM , P I VOT . I NDEX ♦ I P I VDT ) 

C 

DO  11  I»1.3 

11  CONST(  n*FVEC(  I  .1  ) 

C 

DC  20  1=1 *3 
20  XLAM( I  )»C0NST(  I  ) 

PRINT903»XTAU» (XLAM( I) t 1*1 *3) 

903  FORMATdHO/ 

1  IXllHTHICKNESS  =♦  E16.8  / 

2  IXllHALBEDO(X)  =♦  20HA  B*TANH  ( 1 0*  (  X-C  )  )  // 

3  1X3HA  »♦  E16.e»  10X3HB  =♦  E16.8»  10X3HC  «»  E16.8  //) 
C 

CNEW  APPROXIMATION 

c 

M*1 

L»0 

DO  12  I*1»N 
DO  12  J»ltl 
L*L  +  1 

SUM*P(L»M) 

DO  13  K*l,3 

13  SUM  -SUM  4-  CONST<K)#H(L»<»M) 

12  R(L»M)*SUM 
L  =  0 

DO  14  I*1»N 
DO  14  J-1,I 
L'L  +  l 

14  R2( I »J)*R(L»M) 

SIG=0.0 

CALL  OUTPUT 
C 

DO  50  M1»1,M1MAX 
DO  ’  '  M2*l .NPRNT 
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M  =  M+  1 
1  =  0 

DO  15  I = 1 . N 
DO  15  J=1 . r 
L  =  L  +  1 

SUM=P(L.M) 

DO  16  K=l»3 

16  SUM  =  SUM  +  CONr,T(<)#H,L.K;»M) 

15  R(L.M)=SUM 

L  =  C 

DO  17  I=1,N 
DO  17  J=1,I 
L  =  L+1 

17  R2( I »J)=R(L»M) 

18  SIG=-SIG  +  DELTA 
50  CALL  OUTPUT 

C 

RETURN 

END 

SI8FTC  OUTPUT 

SUBROUTINE  OUTPUT 
DIMENSION  XOi 

COMMON  N.RT ( 7 ) .WT ( 7 ) .WR ( 7 ) .AR ( 7.7 ) . NPPNT ,M1 MAX . <MAX . DEL TA . XTAU » 

1  XLAM(3).  B2<7,7) .R2(7.7).IFLAG.R(28.101).T(1491).SIG. 

2  P(28.101)»H(28.3.101) .PLAM(3) .HLAM(3,3)  ♦P2(7»7)  . 

3  H2( 7,7,3) jCONSTI 3) .NEQ 
DO  1  1=1, N 

DO  1  J-.I,N 
1  R2 ( I » J)»R2( J,  I ) 

Y*XTAU*SIG 
DC  3  1=1,3 
3  X(I)»XLAM(I) 

CALL  ALBE00( Y,x,Z) 

PRINTIOC,  SIG,Y,Z 

100  FORMATdHO  7HSIGMA  =,F6.2»  4X5HTAU  =,  F6,2,  4X8HALBEOO  «,F6.2/) 
DO  2  J»1,N 

2  PRINTl0l,J,(R2( I,J) »I=l,N) 

101  FORMATdlO,  7F10.£.) 

RETURN 

END 

sibftc  albedo 

SUBROUTINE  ALBEDO ( Y , X , Z  ) 

DIMENSION  X(3) 

COMMON  N.RT ( 7) ,WT( 7 ) ,WR( 7) ,AR ( 7,7 ) , NPRNT , Ml  MAX , KMA X , DELTA , XTAU , 

1  XLAM{3),  B2  (7,7) ,R2(7,7) , IFLAG,R{28,101 ) ,T( 1491 ) ,SIG, 

2  “P(28,101),H(28,3,101) ,PLAM{3) ,HLAM.(3,3) ,P2(7,7) , 

3  H2(7,7,3) ,CONST(3)»NEQ 
ARG=iO.0*(Y-X{3)  ) 

Z  =  Xd)  +  X(2)*TANH{ARG) 

RETURN 

END 
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